TASR Analysis
Gordon Gianniny
2023-06-26
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6 ✔ purrr 0.3.4
## ✔ tibble 3.2.1 ✔ dplyr 1.1.2
## ✔ tidyr 1.2.0 ✔ stringr 1.4.1
## ✔ readr 2.1.2 ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
##
## Attaching package: 'lubridate'
##
##
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
##
##
## corrplot 0.92 loaded
##
##
## Attaching package: 'ggpp'
##
##
## The following object is masked from 'package:ggplot2':
##
## annotateOverview
This script reads in cleaned datafiles with stream temperature and invertebrate data produced by the invert_data_cleaning.Rmd and teton_data_cleaning.Rmd scripts and performs more involved analyses of stream temperature and invertebrate density/diversity over time. A few findings of interest from previous analyses are that:
* Average summer (July-Aug) stream temperatures have not changed significantly over time
* Average summer (July-Aug) air temperatures have increased over time
* Overall, invetebrate richness has not changed over time
* Invertebrate communities in snow-fed streams appear to be more variable than those in glacier-fed streamsTo build on these findings, the main sections of this script are:
I. Temperature:
1. Overall trends in stream temperature over time - Time series analysis: Site-by-site and grouped by source
2. Trends in summer mean, max, min, and degree-days over time: Site-by-site and grouped by source
4. Modeling stream temperature using SNOTEL data: Site-by-site and grouped by sourceII. Invertebrates:
1. Trends in richness, abundance, and density over time: Site-by-site and grouped by source
2. Modeling richness, abundance, and density using SNOTEL data: Site-by-site and grouped by sourceIII. Temperature and Invertebrates
1. Test for correlations between temperature metrics (mean/max/min summer; degree days): Site-by-site and grouped by sourceRead in necessary files:
stream_temp <- read.csv("Temperature/cleaned_full_datasets/temps_hourly.csv") #hourly stream temp data
invert_density <- read.csv("invert_data/cleaned_csv/full_invert_densities.csv")%>%
rename(site = Stream)#invert density data
invert_richness <- read.csv("invert_data/cleaned_csv/invert_richness.csv") #invert richness data
snotel <- read.csv("teton_snotel.csv") #snotel data
sources <- read.csv("source_info.csv")%>%
rename(site = stream) #source info, rename for merge
### adding source info to stream temps and invert datasets:
stream_source <- merge(stream_temp, sources, all = T)%>%
mutate(date1 = ymd(date1)) #Add source info to hourly temperature data; re-code date as r date format
stream_source_daily <- stream_source%>%
group_by(site, full_name, stream_code, date1, year)%>%
summarise(t_xbar = mean(temp_c, na.rm = T),
t_max = max(temp_c, na.rm = T),
t_min = min(temp_c, na.rm = T),
source = unique(source)) #calculate daily mean, min, and max temps for each site## Warning: There were 10690 warnings in `summarise()`.
## The first warning was:
## ℹ In argument: `t_max = max(temp_c, na.rm = T)`.
## ℹ In group 41: `site = "cloudveil"`, `full_name = "Cloudveil"`, `stream_code =
## "CV"`, `date1 = 2019-09-11`, `year = NA`.
## Caused by warning in `max()`:
## ! no non-missing arguments to max; returning -Inf
## ℹ Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 10689 remaining warnings.## `summarise()` has grouped output by 'site', 'full_name', 'stream_code',
## 'date1'. You can override using the `.groups` argument.density_source <- merge(invert_density, sources, all = T)
richness_source <- merge(invert_richness, sources, all = T)I. Temperature:
1. Trends in summer mean, max, min, and degree-days over time
A. Site-by-site
First, calculate summer (july-aug) mean daily temp, avg daily max and min; plus degree days over 15c? for each site:
stream_summer <- stream_source_daily %>%
mutate(mo = month(date1))%>% #add month column
filter(mo == 7 | mo == 8)%>% #extract July and August
group_by(site, full_name, stream_code, year)%>% #group by site and year
summarise(summer_xbar = mean(t_xbar, na.rm = T),
max_xbar = mean(t_max, na.rm = T),
min_xbar = mean(t_min, na.rm= T), #calculate mean temperature, mean daily max and min
dd_2.5 = length(which(t_max >= 2.5)),
dd_5 = length(which(t_max >= 5)),
dd_10 = length(which(t_max >= 10)), #calculate no. days with max temp >2.5, 5, and 10 deg
source = unique(source)) #retain source col## `summarise()` has grouped output by 'site', 'full_name', 'stream_code'. You can
## override using the `.groups` argument.Next, extract sites with at least 3 years of data:
under3 <- stream_summer %>%
group_by(site)%>% #group by site
filter(!is.na(year))%>% #remove rows with no year
summarise(n_yrs = length(unique(year)))%>% #count num unique years for site
filter(n_yrs < 3) #extract sites with <3yrs data
summer_3plus <- stream_summer %>%
filter(!site%in%under3$site)Trends in mean temp:
Nest data by site:
summer_nested <- summer_3plus %>%
filter(!is.na(year))%>%
nest(data = -site)lm for each site:
summer.mean.lms <- summer_nested %>%
mutate(model = map(data, ~lm(summer_xbar~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
summer.mean.lms## # A tibble: 11 × 7
## # Groups: site [11]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 cloudveil <gropd_df [3 × 10]> year 0.458 0.309 1.48 0.378
## 2 delta <gropd_df [4 × 10]> year 0.0246 0.0193 1.28 0.330
## 3 grizzly <gropd_df [4 × 10]> year 1.21 0.716 1.68 0.234
## 4 mid_teton <gropd_df [6 × 10]> year 0.0296 0.0280 1.05 0.351
## 5 n_teton <gropd_df [6 × 10]> year 0.255 0.460 0.554 0.609
## 6 paintbrush <gropd_df [5 × 10]> year 1.06 0.510 2.08 0.129
## 7 s_ak_basin <gropd_df [4 × 10]> year 0.0539 0.0159 3.40 0.0768
## 8 s_cascade <gropd_df [3 × 10]> year -0.0705 0.0549 -1.28 0.421
## 9 s_teton <gropd_df [5 × 10]> year 0.845 0.422 2.00 0.139
## 10 skillet <gropd_df [3 × 10]> year 0.370 0.0617 5.99 0.105
## 11 windcave <gropd_df [3 × 10]> year 0.119 0.0672 1.77 0.328Nothing significant at 0.05 level; weakly significant increase in mean summer temp at S AK basin (p = 0.08)
Visualization:
summer_mean_withps <- merge(summer_3plus, summer.mean.lms)
pres.pal <- c("#C26ED6", "#F89225", "#76D96F")
summer.mean <- ggplot(summer_mean_withps, aes(x = year, y = summer_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.1, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "Mean daily July-Aug. stream temp., C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
summer.mean## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 11 rows containing non-finite values (stat_smooth).## Warning: Removed 11 rows containing missing values (geom_point).## Warning: Removed 3 row(s) containing missing values (geom_path).
Mean daily max temp
summer.max.lms <- summer_nested %>%
mutate(model = map(data, ~lm(max_xbar~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
summer.max.lms## # A tibble: 11 × 7
## # Groups: site [11]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 cloudveil <gropd_df [3 × 10]> year 0.532 0.328 1.62 0.352
## 2 delta <gropd_df [4 × 10]> year 0.0654 0.0816 0.801 0.507
## 3 grizzly <gropd_df [4 × 10]> year 1.66 1.10 1.51 0.271
## 4 mid_teton <gropd_df [6 × 10]> year -0.0257 0.0485 -0.529 0.625
## 5 n_teton <gropd_df [6 × 10]> year 0.228 0.603 0.378 0.724
## 6 paintbrush <gropd_df [5 × 10]> year 1.74 0.853 2.04 0.135
## 7 s_ak_basin <gropd_df [4 × 10]> year 0.0774 0.108 0.713 0.550
## 8 s_cascade <gropd_df [3 × 10]> year -0.130 0.133 -0.981 0.506
## 9 s_teton <gropd_df [5 × 10]> year 1.16 0.636 1.82 0.166
## 10 skillet <gropd_df [3 × 10]> year 0.365 0.0254 14.4 0.0442
## 11 windcave <gropd_df [3 × 10]> year 0.158 0.101 1.56 0.364Significant increase at Skillet, non-significant increases at all but S cascade and mid teton
Visualization:
summer_max_withps <- merge(summer_3plus, summer.max.lms)
summer.max <- ggplot(summer_max_withps, aes(x = year, y = max_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.1, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "Mean daily max. July-Aug. stream temp., C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
summer.max## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 11 rows containing non-finite values (stat_smooth).## Warning: Removed 11 rows containing missing values (geom_point).## Warning: Removed 3 row(s) containing missing values (geom_path).
Mean daily min temp
summer.min.lms <- summer_nested %>%
mutate(model = map(data, ~lm(min_xbar~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
summer.min.lms## # A tibble: 11 × 7
## # Groups: site [11]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 cloudveil <gropd_df [3 × 10]> year 0.400 0.273 1.46 0.382
## 2 delta <gropd_df [4 × 10]> year 0.0156 0.0152 1.02 0.413
## 3 grizzly <gropd_df [4 × 10]> year 0.849 0.493 1.72 0.227
## 4 mid_teton <gropd_df [6 × 10]> year 0.0639 0.0261 2.45 0.0708
## 5 n_teton <gropd_df [6 × 10]> year 0.254 0.334 0.761 0.489
## 6 paintbrush <gropd_df [5 × 10]> year 0.662 0.325 2.04 0.135
## 7 s_ak_basin <gropd_df [4 × 10]> year 0.0398 0.0440 0.906 0.461
## 8 s_cascade <gropd_df [3 × 10]> year -0.0249 0.0117 -2.14 0.279
## 9 s_teton <gropd_df [5 × 10]> year 0.573 0.290 1.98 0.143
## 10 skillet <gropd_df [3 × 10]> year 0.347 0.0963 3.61 0.172
## 11 windcave <gropd_df [3 × 10]> year 0.0993 0.0354 2.80 0.218Non-significant increases at all sites except south cascade
Visualization:
summer_min_withps <- merge(summer_3plus, summer.min.lms)
summer.min <- ggplot(summer_min_withps, aes(x = year, y = min_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.1, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "Mean daily min. July-Aug. stream temp., C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
summer.min## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 11 rows containing non-finite values (stat_smooth).## Warning: Removed 11 rows containing missing values (geom_point).## Warning: Removed 3 row(s) containing missing values (geom_path).
Degree Days
Days over 2.5 C:
summer.2.5.lms <- summer_nested %>%
mutate(model = map(data, ~lm(dd_2.5~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `model = map(data, ~lm(dd_2.5 ~ year, data = .) %>% tidy)`.
## ℹ In group 8: `site = "s_cascade"`.
## Caused by warning in `summary.lm()`:
## ! essentially perfect fit: summary may be unreliablesummer.2.5.lms## # A tibble: 11 × 7
## # Groups: site [11]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 cloudveil <gropd_df [3 × 10]> year 6.00e+ 0 6.35e+ 0 0.945 0.518
## 2 delta <gropd_df [4 × 10]> year 4.43e+ 0 2.90e+ 0 1.53 0.266
## 3 grizzly <gropd_df [4 × 10]> year 6.1 e+ 0 5.16e+ 0 1.18 0.359
## 4 mid_teton <gropd_df [6 × 10]> year 1.21e+ 0 2.13e+ 0 0.569 0.600
## 5 n_teton <gropd_df [6 × 10]> year 1.90e+ 0 2.58e+ 0 0.736 0.503
## 6 paintbrush <gropd_df [5 × 10]> year 5.70e+ 0 3.70e+ 0 1.54 0.221
## 7 s_ak_basin <gropd_df [4 × 10]> year 5.71e- 1 2.05e+ 0 0.279 0.806
## 8 s_cascade <gropd_df [3 × 10]> year -1.09e-15 1.72e-16 -6.35 0.0994
## 9 s_teton <gropd_df [5 × 10]> year 4.80e+ 0 2.17e+ 0 2.22 0.114
## 10 skillet <gropd_df [3 × 10]> year 5.79e+ 0 2.85e+ 0 2.03 0.291
## 11 windcave <gropd_df [3 × 10]> year 1.55e+ 0 8.25e+ 0 0.188 0.882Non-significant increases at all sites except south cascade
Visualization:
summer_2.5_withps <- merge(summer_3plus, summer.2.5.lms)
summer.2.5 <- ggplot(summer_2.5_withps, aes(x = year, y = dd_2.5, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.1, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "No. summer (July-Aug.) days with max temp >2.5 C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
summer.2.5## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 11 rows containing non-finite values (stat_smooth).## Warning: Removed 11 rows containing missing values (geom_point).## Warning: Removed 3 row(s) containing missing values (geom_path).
Days over 5 C:
summer.5.lms <- summer_nested %>%
mutate(model = map(data, ~lm(dd_5~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
summer.5.lms## # A tibble: 11 × 7
## # Groups: site [11]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 cloudveil <gropd_df [3 × 10]> year 2.50 1.44 1.73 0.333
## 2 delta <gropd_df [4 × 10]> year 0 0 NaN NaN
## 3 grizzly <gropd_df [4 × 10]> year 6.7 5.08 1.32 0.318
## 4 mid_teton <gropd_df [6 × 10]> year -0.0286 0.188 -0.152 0.887
## 5 n_teton <gropd_df [6 × 10]> year 2.09 2.33 0.896 0.421
## 6 paintbrush <gropd_df [5 × 10]> year 9.30 4.66 2.00 0.140
## 7 s_ak_basin <gropd_df [4 × 10]> year 0.357 1.27 0.282 0.805
## 8 s_cascade <gropd_df [3 × 10]> year -1.61 2.46 -0.656 0.630
## 9 s_teton <gropd_df [5 × 10]> year 5.10 2.40 2.12 0.124
## 10 skillet <gropd_df [3 × 10]> year 6.21 2.35 2.64 0.230
## 11 windcave <gropd_df [3 × 10]> year 1.68 0.729 2.31 0.260Non-significant increases at all sites except south cascade, mid teton, and delta (0 all years)
Visualization:
summer_5_withps <- merge(summer_3plus, summer.5.lms)
summer.5 <- ggplot(summer_5_withps, aes(x = year, y = dd_5, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.1, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "No. summer (July-Aug.) days with max temp >5 C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
summer.5## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 11 rows containing non-finite values (stat_smooth).## Warning: Removed 11 rows containing missing values (geom_point).## Warning: Removed 3 row(s) containing missing values (geom_path).
Days over 10 C:
summer.10.lms <- summer_nested %>%
mutate(model = map(data, ~lm(dd_10~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
summer.10.lms## # A tibble: 11 × 7
## # Groups: site [11]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 cloudveil <gropd_df [3 × 10]> year 0 0 NaN NaN
## 2 delta <gropd_df [4 × 10]> year 0 0 NaN NaN
## 3 grizzly <gropd_df [4 × 10]> year 6.1 e+ 0 6.01 1.02e+ 0 0.417
## 4 mid_teton <gropd_df [6 × 10]> year 0 0 NaN NaN
## 5 n_teton <gropd_df [6 × 10]> year 2.01e+ 0 3.20 6.30e- 1 0.563
## 6 paintbrush <gropd_df [5 × 10]> year 5.30e+ 0 3.70 1.43e+ 0 0.247
## 7 s_ak_basin <gropd_df [4 × 10]> year -5.93e-17 0.327 -1.81e-16 1
## 8 s_cascade <gropd_df [3 × 10]> year 0 0 NaN NaN
## 9 s_teton <gropd_df [5 × 10]> year 5.10e+ 0 3.48 1.47e+ 0 0.239
## 10 skillet <gropd_df [3 × 10]> year 0 0 NaN NaN
## 11 windcave <gropd_df [3 × 10]> year 0 0 NaN NaNNothing significant.
Visualization:
summer_10_withps <- merge(summer_3plus, summer.10.lms)
summer.10 <- ggplot(summer_10_withps, aes(x = year, y = dd_10, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.1, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "No. summer (July-Aug.) days with max temp >10 C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
summer.10## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 11 rows containing non-finite values (stat_smooth).## Warning: Removed 11 rows containing missing values (geom_point).## Warning: Removed 3 row(s) containing missing values (geom_path).
B. Grouped by source
First, group by source and re-calculate max, mean, min, and degree days:
summer_source <- stream_summer %>%
group_by(source, year)%>% # group by source type and year
summarise(t_xbar = mean(summer_xbar),
t_max = mean(max_xbar),
t_min = mean(min_xbar), #calculate mean daily temp, mean daily min and max
dd2.5_xbar = mean(dd_2.5),
dd5_xbar = mean(dd_5),
dd10_xbar = mean(dd_10))%>% #calculate mean no. days with max temp >2.5, 5, and 10 C.
filter(!is.na(year), !is.na(source))## `summarise()` has grouped output by 'source'. You can override using the
## `.groups` argument.i. Are there differences in min, max, avg temps between sources?
Visualization:
source.means <- ggplot(summer_source, aes(x = source, y = t_xbar, fill = source))+
geom_boxplot()+
theme_classic()+
scale_fill_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Subterranean Ice"))+
theme(legend.position = "none")+
labs(x = "", y = "Mean daily summer (July-Aug) temp., C")
source.max <- ggplot(summer_source, aes(x = source, y = t_max, fill = source))+
geom_boxplot()+
theme_classic()+
scale_fill_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Subterranean Ice"))+
theme(legend.position = "none")+
labs(x = "Source", y = "Mean daily max. summer (July-Aug) temp., C")
source.min <- ggplot(summer_source, aes(x = source, y = t_max, fill = source))+
geom_boxplot()+
theme_classic()+
scale_fill_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Subterranean Ice"))+
theme(legend.position = "none")+
labs(x = "", y = "Mean daily min. summer (July-Aug) temp., C")
source.means|source.max|source.min
Snow melt streams look to be significantly warmer and have higher max and mins - no surprise there. Check significance of these differences:
max.aov <- aov(t_xbar~source, data = summer_source)
summary(max.aov)## Df Sum Sq Mean Sq F value Pr(>F)
## source 2 55.39 27.693 23.29 2.51e-05 ***
## Residuals 15 17.84 1.189
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(max.aov)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = t_xbar ~ source, data = summer_source)
##
## $source
## diff lwr upr p adj
## snowmelt-glacier 3.67758521 2.042235 5.312935 0.0000909
## sub_ice-glacier -0.08553797 -1.720888 1.549812 0.9898821
## sub_ice-snowmelt -3.76312318 -5.398473 -2.127773 0.0000711Snowmelt significantly warmer than sub-ice or glacier, no difference between sub-ice and glacier fed.
ii. LMs of mean/min/max over time:
First, visualization - re-organize data for easier plotting:
source_plot <- summer_source %>%
select(source, year, t_xbar, t_max, t_min)%>%
pivot_longer(-c(source, year), names_to = "measure", values_to = "temp")Plotting:
source.summer <- ggplot(source_plot, aes(x = year, y = temp, color = measure))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~source, scales = "free_x")+
scale_color_manual(values = c("darkred","darkblue", "darkgreen"), labels = c("Max. temp.", "Min. temp.", "Mean temp"))+
labs(x = "Date", y = "Stream temp, C.", color = "Measurement Type")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
source.summer## `geom_smooth()` using formula 'y ~ x'
Test for significance:
Nest plot data by source and measurement type:
source_nested <- source_plot%>%
nest(data = -c(source, measure))Modeling:
source.lms <- source_nested %>%
mutate(model = map(data, ~lm(temp~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
source.lms## # A tibble: 9 × 8
## # Groups: source [3]
## source measure data term estimate std.error statistic p.value
## <chr> <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 glacier t_xbar <tibble [6 × 2]> year 0.175 0.0225 7.77 0.00148
## 2 glacier t_max <tibble [6 × 2]> year 0.134 0.0514 2.61 0.0594
## 3 glacier t_min <tibble [6 × 2]> year 0.181 0.0283 6.38 0.00309
## 4 snowmelt t_xbar <tibble [6 × 2]> year 0.130 0.408 0.319 0.766
## 5 snowmelt t_max <tibble [6 × 2]> year 0.314 0.531 0.592 0.586
## 6 snowmelt t_min <tibble [6 × 2]> year 0.0448 0.313 0.143 0.893
## 7 sub_ice t_xbar <tibble [6 × 2]> year -0.0517 0.0986 -0.524 0.628
## 8 sub_ice t_max <tibble [6 × 2]> year -0.000823 0.131 -0.00629 0.995
## 9 sub_ice t_min <tibble [6 × 2]> year -0.0694 0.107 -0.647 0.553Significant increases in min, mean, and max summer temp in glacier-fed streams, no change in snowmelt fed streams, non-significant increases in sub_ice streams.
Adding these P-values to the plot:
source_pvals <- merge(source.lms, summer_source)%>% #add model info to source dataframe
mutate(ypos = case_when(measure == "t_max" ~ 0.95, measure == "t_xbar" ~ 0.9, measure == "t_min" ~ 0.85)) #create column for adjusting P-value label Y position
source.summer.pv<- source.summer+
geom_text_npc(data = source_pvals, aes(npcx = 0.1, npcy = ypos, color = measure, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
source.summer.pv## `geom_smooth()` using formula 'y ~ x'
Save:
ggsave("Temperature/plots/source_summer.pdf", source.summer.pv, device = "pdf", units = "in", height = 5, width = 6.5, dpi = "retina")## `geom_smooth()` using formula 'y ~ x'iii. LMs of degree days over time:
Re-organize data for easier plotting and modeling:
sourcedd_plot <- summer_source %>%
select(source, year, dd2.5_xbar, dd5_xbar, dd10_xbar)%>%
pivot_longer(-c(source, year), names_to = "dd_type", values_to = "no_days")Test for significance:
Nest plot data by source and measurement type:
sourcedd_nested <- sourcedd_plot%>%
nest(data = -c(source, dd_type))Modeling:
sourcedd.lms <- sourcedd_nested %>%
mutate(model = map(data, ~lm(no_days~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
sourcedd.lms## # A tibble: 9 × 8
## # Groups: source [3]
## source dd_type data term estimate std.error statistic p.value
## <chr> <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 glacier dd2.5_xbar <tibble> year 2.43 1.30 1.86 0.136
## 2 glacier dd5_xbar <tibble> year 1.79 1.30 1.38 0.241
## 3 glacier dd10_xbar <tibble> year 0 0 NaN NaN
## 4 snowmelt dd2.5_xbar <tibble> year -2.11 4.06 -0.519 0.631
## 5 snowmelt dd5_xbar <tibble> year -1.28 3.93 -0.325 0.762
## 6 snowmelt dd10_xbar <tibble> year -0.671 3.30 -0.204 0.849
## 7 sub_ice dd2.5_xbar <tibble> year -0.518 3.32 -0.156 0.883
## 8 sub_ice dd5_xbar <tibble> year 1.00 1.04 0.963 0.390
## 9 sub_ice dd10_xbar <tibble> year 0.0714 0.169 0.423 0.694No significant changes in degree days in any group.
Plotting with these P-values:
sourcedd_pvals <- merge(sourcedd.lms, sourcedd_plot)%>% #add model info to source dataframe
mutate(ypos = case_when(dd_type == "dd2.5_xbar" ~ 0.85, dd_type == "dd5_xbar" ~ 0.9, dd_type == "dd10_xbar" ~ 0.95)) #create column for adjusting P-value label Y position
sourcedd.summer<- ggplot(sourcedd_pvals, aes(x = year, y = no_days, color = dd_type))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~source, scales = "free_x")+
scale_color_manual(values = c("darkred","orange", "darkorange"), labels = c(">10 C", ">2.5 C", ">5 C"))+
labs(x = "Year", y = "Num. degree days", color = "Degree day threshold")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))+
geom_text_npc(aes(npcx = 0.9, npcy = ypos, color = dd_type, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
sourcedd.summer## `geom_smooth()` using formula 'y ~ x'
iv. Lumping glacier & snowmelt:
First, combine glacier and snowmelt streams and recalculate mean, min, and max summer temperatures:
lumped_source <- stream_summer %>%
mutate(source_lumped = ifelse(source == "sub_ice", "sub_ice", "other"))%>% #add new col with lumped sources
filter(!is.na(year), !is.na(source))%>% #remove rows with no year or source
group_by(source_lumped, year)%>% #group by lumped source and year
summarise(t_xbar = mean(summer_xbar), #calculate mean, min, and max summer temps
max_xbar = mean(max_xbar),
min_xbar = mean(min_xbar))%>%
pivot_longer(!c(source_lumped, year), names_to = "measure", values_to = "temp") #pivot longer for plotting + analysis## `summarise()` has grouped output by 'source_lumped'. You can override using the
## `.groups` argument.Nest data by lumped source and measurement type:
lumped_nested <- lumped_source%>%
nest(data = -c(source_lumped, measure))lm for each source + year combo:
lumped.lms <- lumped_nested %>%
mutate(model = map(data, ~lm(temp~year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "year") #remove intercept term to just get info on year for each site
lumped.lms## # A tibble: 6 × 8
## # Groups: source_lumped [2]
## source_lumped measure data term estimate std.error statistic p.value
## <chr> <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 other t_xbar <tibble> year 0.196 0.264 0.745 0.498
## 2 other max_xbar <tibble> year 0.295 0.356 0.829 0.453
## 3 other min_xbar <tibble> year 0.141 0.198 0.709 0.518
## 4 sub_ice t_xbar <tibble> year -0.0517 0.0986 -0.524 0.628
## 5 sub_ice max_xbar <tibble> year -0.000823 0.131 -0.00629 0.995
## 6 sub_ice min_xbar <tibble> year -0.0694 0.107 -0.647 0.553Add p-values to data and plot:
lumped_pvals <- merge(lumped_source, lumped.lms)%>% #add model info to source dataframe
mutate(ypos = case_when(measure == "t_xbar" ~ 0.9, measure == "max_xbar" ~ 0.95, measure == "min_xbar" ~ 0.85)) #create column for adjusting P-value label Y position
lumped.summer<- ggplot(lumped_pvals, aes(x = year, y = temp, color = measure))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~source_lumped, scales = "free_x")+
scale_color_manual(values = c("darkred","darkblue", "darkgreen"), labels = c("Max. temp.", "Min. temp.", "Mean temp"))+
labs(x = "Year", y = "Stream temperature, C", color = "Measurement Type")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))+
geom_text_npc(aes(npcx = 0.1, npcy = ypos, color = measure, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
lumped.summer## `geom_smooth()` using formula 'y ~ x'
ggsave("Temperature/plots/lumped_summer.pdf", lumped.summer, device = "pdf", units = "in", height = 4, width = 6.5, dpi = "retina")## `geom_smooth()` using formula 'y ~ x'2. Modeling stream temperature using SNOTEL data
First, calculate April 1 SWE and summer air temp from snotel data; add to mean temperature dataset:
#Calc means across both sites:
snotel_means <- snotel %>%
mutate(date1 = as.Date(date, format = "%m/%d/%y"), #recode date column as R date
airtemp_c = as.numeric(airtemp_c))%>% #recode airtemp as numeric
group_by(date1)%>% #group by date
summarise(swe_xbar = mean(swe_cm), #calculate mean SWE and airtemp across both stations
airtemp_xbar = mean(airtemp_c, na.rm = T))%>%
mutate(mo = month(date1), yr = year(date1)) #Add month and year cols## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `airtemp_c = as.numeric(airtemp_c)`.
## Caused by warning:
## ! NAs introduced by coercion#Calculate max SWE and mean summer airtemp for each year:
snotel_summary <- snotel_means %>%
group_by(yr)%>% #group by year
mutate(swe_max = max(swe_xbar))%>% #extract max SWE and add to new col
ungroup()%>% #ungroup to get all columns back
filter(mo == 7 | mo == 8)%>% #extract summer months
group_by(yr)%>% #group by year
summarise(airtemp_summer = mean(airtemp_xbar), #calculate mean summer air temp
swe_max = unique(swe_max))%>% #keep max SWE
rename(year = yr)
#Merge with long-term temperature averages:
stream_snotel <- merge(summer_3plus, snotel_summary)A. Site-by-site
i. Correlations between SWE and Summer Temp:
LMs - first, nest by site:
stream_snotel_nested <- stream_snotel %>%
nest(data = - site)lm for each site:
swe.temp.lms <- stream_snotel_nested %>%
mutate(model = map(data, ~lm(summer_xbar~swe_max, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "swe_max") #remove intercept term to just get info on year for each site
swe.temp.lms## # A tibble: 11 × 7
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 mid_teton <tibble [6 × 12]> swe_max 0.000784 0.00325 0.241 0.821
## 2 windcave <tibble [3 × 12]> swe_max -0.00146 0.0127 -0.115 0.927
## 3 n_teton <tibble [6 × 12]> swe_max -0.0561 0.0405 -1.39 0.238
## 4 s_cascade <tibble [3 × 12]> swe_max -0.0165 0.00937 -1.77 0.328
## 5 s_ak_basin <tibble [4 × 12]> swe_max 0.00414 0.00603 0.687 0.563
## 6 s_teton <tibble [5 × 12]> swe_max -0.0487 0.0449 -1.08 0.357
## 7 delta <tibble [4 × 12]> swe_max -0.00178 0.00155 -1.15 0.370
## 8 paintbrush <tibble [5 × 12]> swe_max -0.0634 0.0543 -1.17 0.327
## 9 grizzly <tibble [4 × 12]> swe_max -0.0567 0.0840 -0.675 0.569
## 10 skillet <tibble [3 × 12]> swe_max -0.0305 0.00316 -9.67 0.0656
## 11 cloudveil <tibble [3 × 12]> swe_max -0.0166 0.0329 -0.504 0.703Add P values and merge:
swe_summer_ps <- merge(swe.temp.lms, stream_snotel)
swe.summer.temp <- ggplot(swe_summer_ps, aes(x = swe_max, y = summer_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Annual maxium SWE, cm", y = "Mean daily July-Aug. stream temp., C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
swe.summer.temp## `geom_smooth()` using formula 'y ~ x'
ii. Summer stream temp ~ summer air temp:
lm for each site:
air.temp.lms <- stream_snotel_nested %>%
mutate(model = map(data, ~lm(summer_xbar~airtemp_summer, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "airtemp_summer") #remove intercept term to just get info on year for each site
air.temp.lms## # A tibble: 11 × 7
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 mid_teton <tibble [6 × 12]> airtemp_su… 0.0643 0.0829 0.775 0.481
## 2 windcave <tibble [3 × 12]> airtemp_su… 0.112 0.273 0.410 0.752
## 3 n_teton <tibble [6 × 12]> airtemp_su… -0.457 1.32 -0.346 0.747
## 4 s_cascade <tibble [3 × 12]> airtemp_su… -0.251 0.0964 -2.60 0.233
## 5 s_ak_basin <tibble [4 × 12]> airtemp_su… 0.348 0.250 1.39 0.299
## 6 s_teton <tibble [5 × 12]> airtemp_su… 0.0288 4.21 0.00683 0.995
## 7 delta <tibble [4 × 12]> airtemp_su… -0.214 0.0826 -2.59 0.122
## 8 paintbrush <tibble [5 × 12]> airtemp_su… -0.0494 5.21 -0.00947 0.993
## 9 grizzly <tibble [4 × 12]> airtemp_su… 3.18 5.68 0.559 0.632
## 10 skillet <tibble [3 × 12]> airtemp_su… -0.346 2.58 -0.134 0.915
## 11 cloudveil <tibble [3 × 12]> airtemp_su… 1.81 1.04 1.74 0.332Add P values and plot:
air_summer_ps <- merge(air.temp.lms, stream_snotel)
swe.summer.temp <- ggplot(air_summer_ps, aes(x = airtemp_summer, y = summer_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Mean daily July-Aug. air temp, C", y = "Mean daily July-Aug. stream temp., C", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
swe.summer.temp## `geom_smooth()` using formula 'y ~ x'
B. Grouped by source
First, add SNOTEL summary info to mean summer temps grouped by source:
source_snotel <- merge(summer_source, snotel_summary)i. LMs of mean/min/max temp ~ SWE:
Re-organize data for analysis and plotting:
source_swe <- source_snotel %>%
select(source, swe_max, t_xbar, t_min, t_max)%>% #select temperature, source, and swe cols
pivot_longer(!c(source, swe_max), names_to = "measure", values_to = "temp") #pivot longer for plotting + analysisNest data by source and measure:
sourceSwe_nested <- source_swe %>%
nest(data = -c(source, measure))LMs for each source and measure:
source.swe.lms <- sourceSwe_nested %>%
mutate(model = map(data, ~lm(temp~swe_max, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "swe_max") #remove intercept term to just get info on year for each site
source.swe.lms## # A tibble: 9 × 8
## source measure data term estimate std.error statistic p.value
## <chr> <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 glacier t_xbar <tibble [6 × 2]> swe_max -0.00523 0.00895 -0.585 0.590
## 2 glacier t_min <tibble [6 × 2]> swe_max -0.00349 0.00961 -0.363 0.735
## 3 glacier t_max <tibble [6 × 2]> swe_max -0.00689 0.00801 -0.861 0.438
## 4 snowmelt t_xbar <tibble [6 × 2]> swe_max -0.0491 0.0349 -1.41 0.231
## 5 snowmelt t_min <tibble [6 × 2]> swe_max -0.0415 0.0249 -1.67 0.170
## 6 snowmelt t_max <tibble [6 × 2]> swe_max -0.0548 0.0502 -1.09 0.336
## 7 sub_ice t_xbar <tibble [6 × 2]> swe_max 0.0192 0.00656 2.92 0.0431
## 8 sub_ice t_min <tibble [6 × 2]> swe_max 0.0231 0.00566 4.07 0.0152
## 9 sub_ice t_max <tibble [6 × 2]> swe_max 0.0106 0.0139 0.765 0.487Add p-values and plot:
sourceSwe_pval<-merge(source_swe, source.swe.lms)%>%
mutate(ypos = case_when(measure == "t_xbar" ~ 0.9, measure == "t_max" ~ 0.95, measure == "t_min" ~ 0.85)) #create column for adjusting P-value label Y position
source.swe<- ggplot(sourceSwe_pval, aes(x = swe_max, y = temp, color = measure))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~source, scales = "free_x")+
scale_color_manual(values = c("darkred","darkblue", "darkgreen"), labels = c("Max. temp.", "Min. temp.", "Mean temp"))+
labs(x = "Annual maximum SWE, cm", y = "Stream temperature, C", color = "Measurement Type")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))+
geom_text_npc(aes(npcx = 0.9, npcy = ypos, color = measure, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
source.swe## `geom_smooth()` using formula 'y ~ x'
Save:
ggsave("Temperature/plots/source_swe.pdf", source.swe, device = "pdf", units = "in", height = 4, width = 6.5, dpi = "retina")## `geom_smooth()` using formula 'y ~ x'ii. LMs of mean/min/max temp ~ Air temp:
Re-organize data for analysis and plotting:
source_air <- source_snotel %>%
select(source, airtemp_summer, t_xbar, t_min, t_max)%>% #select temperature, source, and swe cols
pivot_longer(!c(source, airtemp_summer), names_to = "measure", values_to = "temp") #pivot longer for plotting + analysisNest data by source and measure:
sourceAir_nested <- source_air %>%
nest(data = -c(source, measure))LMs for each source and measure:
source.air.lms <- sourceAir_nested %>%
mutate(model = map(data, ~lm(temp~airtemp_summer, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "airtemp_summer") #remove intercept term to just get info on year for each site
source.air.lms## # A tibble: 9 × 8
## source measure data term estimate std.error statistic p.value
## <chr> <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 glacier t_xbar <tibble [6 × 2]> airtem… 0.249 0.221 1.13 0.322
## 2 glacier t_min <tibble [6 × 2]> airtem… 0.300 0.219 1.37 0.243
## 3 glacier t_max <tibble [6 × 2]> airtem… 0.124 0.229 0.542 0.616
## 4 snowmelt t_xbar <tibble [6 × 2]> airtem… -0.701 1.10 -0.635 0.560
## 5 snowmelt t_min <tibble [6 × 2]> airtem… -0.712 0.806 -0.884 0.427
## 6 snowmelt t_max <tibble [6 × 2]> airtem… -0.488 1.53 -0.318 0.766
## 7 sub_ice t_xbar <tibble [6 × 2]> airtem… 0.147 0.315 0.468 0.664
## 8 sub_ice t_min <tibble [6 × 2]> airtem… 0.312 0.321 0.972 0.386
## 9 sub_ice t_max <tibble [6 × 2]> airtem… -0.117 0.410 -0.284 0.790Add p-values and plot:
sourceAir_pval<-merge(source_air, source.air.lms)%>%
mutate(ypos = case_when(measure == "t_xbar" ~ 0.9, measure == "t_max" ~ 0.95, measure == "t_min" ~ 0.85)) #create column for adjusting P-value label Y position
source.air<- ggplot(sourceAir_pval, aes(x = airtemp_summer, y = temp, color = measure))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~source, scales = "free_x")+
scale_color_manual(values = c("darkred","darkblue", "darkgreen"), labels = c("Max. temp.", "Min. temp.", "Mean temp"))+
labs(x = "Mean July-Aug. air temp., C", y = "Stream temperature, C", color = "Measurement Type")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))+
geom_text_npc(aes(npcx = 0.1, npcy = ypos, color = measure, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
source.air## `geom_smooth()` using formula 'y ~ x'
Save:
ggsave("Temperature/plots/source_air.pdf", source.air, device = "pdf", units = "in", height = 4, width = 6.5, dpi = "retina")## `geom_smooth()` using formula 'y ~ x'3. Overall trends in stream temperature over time - Time series analysis
A. Site-by-site
B. Grouped by source
II. Invertebrates:
1. Trends in richness, abundance, and density over time
A. Site-by-site
First, calculate richness, mean density, and total abundance for each site for each year; remove sites with only 1 year of data:
oneYr <- density_source%>%
group_by(site)%>% #group by site
summarise(nYrs = length(unique(Year)))%>% #count no. years
filter(nYrs <3) #extract sites w/ only 1-2 yrs data
invert_means <- density_source %>%
group_by(stream_code, full_name, site, Year)%>% #group by stream and year
summarise(richness = length(unique(taxa_lwr)), #calculate yearly richness, density, and total abundance
density_xbar = mean(density_xbar),
tot_abund = sum(abundance_total),
source = unique(source))%>%
filter(!site%in%oneYr$site) #drop sites with 1-2yrs of data## `summarise()` has grouped output by 'stream_code', 'full_name', 'site'. You can
## override using the `.groups` argument.i. LM’s of Richness ~ year at each site:
First, nest data by site:
invert_nested <- invert_means %>%
nest(data = -site)LM’s of Richness ~ year for each site:
richness.yr.lm <- invert_nested %>%
mutate(model = map(data, ~lm(richness~Year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "Year") #remove intercept term to just get info on year for each site
richness.yr.lm## # A tibble: 11 × 7
## # Groups: site [11]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 s_ak_basin <gropd_df [6 × 7]> Year -1.34e+ 0 0.956 -1.41e+ 0 0.233
## 2 cloudveil <gropd_df [3 × 7]> Year 1.50e+ 0 0.866 1.73e+ 0 0.333
## 3 delta <gropd_df [5 × 7]> Year -1.10e+ 0 0.907 -1.21e+ 0 0.312
## 4 grizzly <gropd_df [5 × 7]> Year 1.00e- 1 0.790 1.27e- 1 0.907
## 5 mid_teton <gropd_df [7 × 7]> Year -8.93e- 1 0.394 -2.27e+ 0 0.0725
## 6 n_teton <gropd_df [7 × 7]> Year -1.54e+ 0 0.886 -1.73e+ 0 0.144
## 7 paintbrush <gropd_df [5 × 7]> Year -8.00e- 1 0.542 -1.48e+ 0 0.236
## 8 s_cascade <gropd_df [6 × 7]> Year -2.93e+ 0 0.852 -3.44e+ 0 0.0264
## 9 skillet <gropd_df [5 × 7]> Year -1.01e-13 1.25 -8.08e-14 1.00
## 10 s_teton <gropd_df [7 × 7]> Year -1.43e+ 0 0.818 -1.75e+ 0 0.141
## 11 windcave <gropd_df [7 × 7]> Year -7.86e- 1 0.340 -2.31e+ 0 0.0686Add p-values and plot:
richness_withps <- merge(invert_means, richness.yr.lm)
pres.pal <- c("#C26ED6", "#F89225", "#76D96F")
richness.yr <- ggplot(richness_withps, aes(x = Year, y = richness, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "Taxonomic richness", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
richness.yr## `geom_smooth()` using formula 'y ~ x'
Significant decline in richness at S Cascade, non-significant declines at most other sites
Save:
ggsave("invert_data/plots/richness_year.pdf", richness.yr, device = "pdf", dpi = "retina", units = "in", height = 5, width = 6.5)## `geom_smooth()` using formula 'y ~ x'ii. LM’s of Density ~ year at each site:
LM’s of Density ~ year for each site:
dens.yr.lm <- invert_nested %>%
mutate(model = map(data, ~lm(density_xbar~Year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "Year") #remove intercept term to just get info on year for each site
dens.yr.lm## # A tibble: 11 × 7
## # Groups: site [11]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 s_ak_basin <gropd_df [6 × 7]> Year 12.9 12.2 1.06 0.351
## 2 cloudveil <gropd_df [3 × 7]> Year 64.6 16.9 3.82 0.163
## 3 delta <gropd_df [5 × 7]> Year 15.9 41.0 0.389 0.724
## 4 grizzly <gropd_df [5 × 7]> Year 6.06 44.9 0.135 0.901
## 5 mid_teton <gropd_df [7 × 7]> Year 42.9 39.5 1.09 0.327
## 6 n_teton <gropd_df [7 × 7]> Year -7.10 18.6 -0.382 0.718
## 7 paintbrush <gropd_df [5 × 7]> Year -12.9 203. -0.0637 0.953
## 8 s_cascade <gropd_df [6 × 7]> Year -2.25 16.2 -0.139 0.896
## 9 skillet <gropd_df [5 × 7]> Year 32.1 30.6 1.05 0.372
## 10 s_teton <gropd_df [7 × 7]> Year 3.68 2.05 1.79 0.133
## 11 windcave <gropd_df [7 × 7]> Year 28.6 19.6 1.46 0.205Add p-values and plot:
density_withps <- merge(invert_means, dens.yr.lm)
pres.pal <- c("#C26ED6", "#F89225", "#76D96F")
dens.yr <- ggplot(density_withps, aes(x = Year, y = density_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "Mean density, individuals/m2", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
dens.yr## `geom_smooth()` using formula 'y ~ x'
No significant trends in average density at any site
Save:
ggsave("invert_data/plots/density_year.pdf", dens.yr, device = "pdf", dpi = "retina", units = "in", height = 5, width = 6.5)## `geom_smooth()` using formula 'y ~ x'iii. LM’s of total abundance ~ year at each site:
LM’s of Density ~ year for each site:
abund.yr.lm <- invert_nested %>%
mutate(model = map(data, ~lm(tot_abund~Year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "Year") #remove intercept term to just get info on year for each site
abund.yr.lm## # A tibble: 11 × 7
## # Groups: site [11]
## site data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 s_ak_basin <gropd_df [6 × 7]> Year 60.5 149. 0.406 0.705
## 2 cloudveil <gropd_df [3 × 7]> Year 333. 133. 2.51 0.242
## 3 delta <gropd_df [5 × 7]> Year -110. 383. -0.287 0.793
## 4 grizzly <gropd_df [5 × 7]> Year 107. 235. 0.455 0.680
## 5 mid_teton <gropd_df [7 × 7]> Year 4.89 190. 0.0258 0.980
## 6 n_teton <gropd_df [7 × 7]> Year -28.2 109. -0.259 0.806
## 7 paintbrush <gropd_df [5 × 7]> Year -383. 347. -1.10 0.350
## 8 s_cascade <gropd_df [6 × 7]> Year -193. 65.0 -2.97 0.0411
## 9 skillet <gropd_df [5 × 7]> Year 269. 281. 0.955 0.410
## 10 s_teton <gropd_df [7 × 7]> Year -137. 96.4 -1.42 0.215
## 11 windcave <gropd_df [7 × 7]> Year 50.6 119. 0.424 0.689Add p-values and plot:
abund_withps <- merge(invert_means, abund.yr.lm)
pres.pal <- c("#C26ED6", "#F89225", "#76D96F")
abund.yr <- ggplot(density_withps, aes(x = Year, y = tot_abund, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~full_name, scales = "free")+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Date", y = "Total abundance", color = "Stream Source")+
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
abund.yr## `geom_smooth()` using formula 'y ~ x'
No significant trends in total abundance at any site.
Save:
ggsave("invert_data/plots/abundance_year.pdf", abund.yr, device = "pdf", dpi = "retina", units = "in", height = 5, width = 6.5)## `geom_smooth()` using formula 'y ~ x'B. Grouped by source
First, group by source and calculate mean richness, density, and abundance:
source_inverts <- invert_means %>%
group_by(source, Year)%>%
summarise(richness_xbar = mean(richness),
density_xbar = mean(density_xbar),
abund_xbar = mean(tot_abund))## `summarise()` has grouped output by 'source'. You can override using the
## `.groups` argument.i. LMs of richness ~ year for each source:
First, nest data by source:
source_inverts_nested <- source_inverts %>%
nest(data = -source)LMs of richness ~ year for each source:
ricness.source.lm <- source_inverts_nested %>%
mutate(model = map(data, ~lm(richness_xbar~Year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "Year") #remove intercept term to just get info on year for each site
ricness.source.lm## # A tibble: 3 × 7
## # Groups: source [3]
## source data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 glacier <tibble [7 × 4]> Year -0.310 0.341 -0.909 0.405
## 2 snowmelt <tibble [7 × 4]> Year -1.69 0.895 -1.89 0.117
## 3 sub_ice <tibble [7 × 4]> Year -1.51 0.537 -2.80 0.0379Add p-values and plot:
richnessSource_pv <- merge(source_inverts, ricness.source.lm)
pres.pal <- c("#C26ED6", "#F89225", "#76D96F")
richness.source <- ggplot(richnessSource_pv, aes(x = Year, y = richness_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~source)+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.95, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Year", y = "Mean taxonomic richness", color = "Stream Source")+
theme(legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
richness.source## `geom_smooth()` using formula 'y ~ x'
Significant decline in richness at subterranean ice sites
Save:
ggsave("invert_data/plots/richness_source.pdf", richness.source, device = "pdf", dpi = "retina", units = "in", height = 5, width = 6.5)## `geom_smooth()` using formula 'y ~ x'ii. LMs of density ~ year for each source:
LMs of density ~ year for each source:
dens.source.lm <- source_inverts_nested %>%
mutate(model = map(data, ~lm(density_xbar~Year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "Year") #remove intercept term to just get info on year for each site
dens.source.lm## # A tibble: 3 × 7
## # Groups: source [3]
## source data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 glacier <tibble [7 × 4]> Year -12.1 27.2 -0.444 0.676
## 2 snowmelt <tibble [7 × 4]> Year 0.528 30.4 0.0174 0.987
## 3 sub_ice <tibble [7 × 4]> Year 10.7 12.2 0.875 0.422Add p-values and plot:
densSource_pv <- merge(source_inverts, dens.source.lm)
pres.pal <- c("#C26ED6", "#F89225", "#76D96F")
dens.source <- ggplot(densSource_pv, aes(x = Year, y = density_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~source)+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.95, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Year", y = "Mean density, individuals/m2", color = "Stream Source")+
theme(legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
dens.source## `geom_smooth()` using formula 'y ~ x'
No significant trends in average density in any source.
Save:
ggsave("invert_data/plots/density_source.pdf", dens.source, device = "pdf", dpi = "retina", units = "in", height = 5, width = 6.5)## `geom_smooth()` using formula 'y ~ x'iii. LMs of abundance ~ year for each source:
LMs of abundance ~ year for each source:
abund.source.lm <- source_inverts_nested %>%
mutate(model = map(data, ~lm(abund_xbar~Year, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
filter(term == "Year") #remove intercept term to just get info on year for each site
abund.source.lm## # A tibble: 3 × 7
## # Groups: source [3]
## source data term estimate std.error statistic p.value
## <chr> <list> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 glacier <tibble [7 × 4]> Year -103. 159. -0.649 0.545
## 2 snowmelt <tibble [7 × 4]> Year -110. 85.0 -1.29 0.254
## 3 sub_ice <tibble [7 × 4]> Year -47.4 101. -0.470 0.658Add p-values and plot:
abundSource_pv <- merge(source_inverts, abund.source.lm)
pres.pal <- c("#C26ED6", "#F89225", "#76D96F")
abund.source <- ggplot(abundSource_pv, aes(x = Year, y = abund_xbar, color = source))+
geom_point()+
geom_line()+
geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
theme_classic()+
facet_wrap(~source)+
scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
geom_text_npc(aes(npcx = 0.9, npcy = 0.95, label = paste( "P = ",round(p.value, 2), sep = "")))+
labs(x = "Year", y = "Mean density, individuals/m2", color = "Stream Source")+
theme(legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.title = element_text(size = 11, face = "bold"),
axis.text = element_text(size = 8),
legend.text = element_text(size = 10),
legend.title = element_text(size = 11, face = "bold"),
strip.text = element_text(size = 10, face = "bold"))
abund.source## `geom_smooth()` using formula 'y ~ x'
No significant trends in average abundance in any source.
Save:
ggsave("invert_data/plots/abundance_source.pdf", abund.source, device = "pdf", dpi = "retina", units = "in", height = 5, width = 6.5)## `geom_smooth()` using formula 'y ~ x'